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Abstract 

We consider the accuracy of the space discretization for time-dependent 
problems when a nonuniform mesh is used. We show that many schemes reduce to 
first-order accuracy while a popular finite volume scheme is even inconsistent 
for general grids. This accuracy is based on physical variables. However, 
when accuracy is measured in computational variables then second-order 
accuracy can be obtained. This is meaningful only if the mesh accurately 
reflects the properties of the solution. In addition we analyze the stability 
properties of some improved accurate schemes and show that they also allow for 
larger time steps when Runge-Kutta type methods are used to advance in time. 
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I. INTRODUCTION 


With the latest class of computers it is now possible to routinely perform 
two-dimensional calculations for both the Euler and the compressible Navier- 
Stokes equations. Three-dimensional Euler calculations about simple shapes 
are feasible if a relatively coarse grid is used. Three-dimensional 
calculations based on the thin layer compressible Navier-Stokes calculations 
are just becoming available. The current trend is to use coordinates 
constructed so that the solid surfaces are coordinate lines while the other 
coordinate directions are close to orthogonal (see however [4]). This 
approach simplifies the boundary conditions at the solid surfaces. Moreover, 
whenever boundary layers exist, this approach allows one to refine the mesh 
across the boundary layer while having a relatively coarse mesh parallel to 
the boundary layer. 

The construction of a general curvilinear grid is still not easily 
accomplished. Difficulties occur when one wishes the grids to vary smoothly 
while also being able to concentrate points in certain regions. Since many 
bodies have cusped regions, the meshes frequently are far from regular in 
these regions. In three dimensions these difficulties are compounded. First, 
on present day machines one is still restricted to relatively coarse grids. 
In addition the shape of the bodies can be more complicated in three 
dimensions. One frequently uses a quasi two-dimensional approach which has 
difficulties when the body no longer appears in some two-dimensional slice. 
The result of all these difficulties is that the meshes that are constructed 
are highly distorted in some regions. These distortions appear as high aspect 
ratios and angles quite far from 90° for quadrilateral elements. Even more 
disturbing is the change of these quantities from a cell to its neighboring 
cells (see, for example, Figure 2). 
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In this paper we consider two aspects of the calculations for compressible 
flow. The first is the accuracy of various schemes when used on grids which 
have distorted meshes. In addition we shall also consider some explicit time 
marching techniques that utilize improved space discretizations. We shall 
also discuss finite volume approaches and show that some of the improvements 
also allow larger time steps. 

The effect of nonuniform grids on accuracy has previously been discussed 
by several authors, e.g., Hoffman [7], Mastin [9], Steger [14], and Vinokur 
[22]. Rai et . al . [12] have performed calculations on discontinuous grids 
using a first-order upwind scheme. Pike [10] has also considered upwind 
schemes on general (one-dimensional) meshes and has examined conservation and 
TVD properties. 


II. SPACE DISCRETIZATION 

We shall consider time dependent problems. This also applies to steady 
state problems if we use pseudo time-dependent methods. Thus we shall march a 
time-dependent equation which need not be consistent with the time-dependent 
physics of the process but has the same steady state equation. 

We consider the conservation equations 

w+f +g + h =0. (2.1) 

t x & y z 

We decouple the time integration from the space discretization. Hence we 

will, for now, keep the time variable continuous and only consider the 

approximation in space. Using a finite volume approach we integrate (2.1) in 

a cell n and assuming that n is independent of time we obtain 
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wdxdydz 


+ / 

n 


(f + g + h„ )dxdydz 
x y 2 


0 


Using the divergence theorem this becomes 


( 2 . 2 ) 


d 

dt 


J wdV + j 

n an 


f.ndS = 0 


(2.3) 


where F = (f, g, h) 1 " and 


n is the unit outward normal. Equivalently, 


d 

dt 


J wdv 

n 


f dv 


n 


+ -J- / f.ndS = 


an 


(2.4) 


Let w be the average of w in a cell. Then we have 

w = y f wdv (2.3a) 

v n 

and 

+ y f f.ndS = 0. (2.5b) 

n 

We stress that (2.5) is exact and that no approximations have been made to 
this point. 

To illustrate the situation we consider the two-dimensional case shown in 


Figure 1 
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In this case (2.5) becomes 


w = — f wdxdy (2.6a) 

v n 

V + f + f + f + f (fdy - gdx) = 0. (2.6b) 

AB BC CD DA 

Consider, for example, the integral 

( fdy - gdx. (2.7) 

BC 

At this stage we replace (2.7) with an integration rule. We would prefer that 
our formulas be at least second-order accurate for uniform meshes and so we 
consider two integration rules. The first is the midpoint rule. Thus we have 
(see Figure 1) 
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/ (fdy - gdx) = Ayf(G) - Axg(G) + 0((Ax) 3 + (Ay) 3 ). (2.8) 
BC 

We have replaced the line integral by values of the fluxes at the point G, 
where G is the midpoint measured in terms of x and y. Ax, Ay are the 
changes in x and y respectively between the points B and C. An 
alternative is to use the trapezoidal rule. We then have 


/ fdy - gdx =y 2 [Ay(f(B) - f(C)) - Ax(g(B) + g(C))] 

BC 

(2.9) 

+ 0( (Ax) 3 + (Ay) 3 ) . 


To proceed further we must know the fluxes at the point G (for (2.8)) 
or B and C (for (2.9)). The fluxes f, g are functions of the dependent 
variable w and sometimes also the independent variables. However, w is 
defined in (2.6a) as a cell averaged quantity and is not known at specific 
points. To overcome this difficulty we shall replace w by w. . where 
w. . is a point value of w at some point within the cell. An appropriate 

1 9 J 

choice is to use the centroid of the cell and then 


w. . = i f wdv + 0(V 3 ) 

1,3 v h 


( 2 . 10 ) 


and we retain second-order accuracy. 

Having defined w- • we can then find w at the point G or at B 

9 J 

and C by averaging. One possibility (e.g., [8]) is to use simple arithmetic 
averaging so that 


w(G) 


= 72 (w . . + w. , . . ) , 

z i,J i+l, J 


(2.11a) 
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Similarly if one uses (2.9) then we have 


w(C) 




+ w 


i+l,j 


+ w 


i»j+l 




(2.11b) 


We note that we have averaged the quantity w. Since f and g are 
nonlinear functions of w we can instead average f and g or else other 
quantities. The disadvantage of averaging f and g is that this decouples 
the even and odd points when (2.11a) is used. The advantage of averaging f 
and g is that in the steady state the fluxes are continuous across shocks if 
the shock is aligned with a mesh line. 

Having made the choice (2.11a) or (2.11b) this completes the space 
discretization of (2.1). Assuming we average the fluxes the choice of (2.11a) 
leads to the two-dimensional scheme 


dw. 


V(ft) 


— — — + V? f (f . . - . + f. . )Ax i, . - (f. . + f. . )Ax. i, 
dt xl i+I, j i,j i+ V2 ,J 1,3 i-l ,J i- V2 ,3 


+ (f . ... + f. . )Ax. i, - (f. . + f. . . )Ax. . i, 
i»J+l i,3 i.J +1 /2 i >J i,J-l i,J~ 72 


+ <g i + I,j + 8 i,j >Ay i+ V 2 ,J - <g i.j + s i-l,j )Ay l-V 2 ,} 


( 2 . 12 ) 


+ + S i,j )Ay i, j+ V2~ <g i.J + e l,j-l )Ay i,J- V 2 1 ‘ ° 


where V(fi) = l / 2 |(y c “ y A )^ x B “ Xq) “ ( x c - x A ^y B “ y D )| ( see Figure 1). 

Assuming that the positions (x, y) are given at the mesh nodes 

(Ax)^ + ^ _. = x(C) - x(B), etc. For a uniform mesh Ax, Ay constant this 


reduces to a centered difference formula 
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(2.13) 


The trapezoidal rule (2.11b) reduces, in the case of a uniform mesh, to 


(f . . . - f . . . ) f . , . - - f , . , 

2 i+l, J i-l.J + i+l, 3-1 i-l, J-l 


S i+l,j+l g i+l,j~ 


,3-1 , o (g i,j + l ~ g i,j~l ) , g i-l,j + l g i-l.,j-l _ n 

2 Ay 2Ay J * 


We shall see later that (2.14) has advantages in terms of the allowable time 
step and in terms of accuracy, while (2.13) has fewer operations and does not 
require any special treatment for the tangential components near boundaries. 

We will discuss the time integration of these systems later in Section 
V. For later purposes we note that there is another possibility besides 

(2.13) and (2.14). This is given by (for uniform meshes) 


dW i,j + 1, £ i+l,j+l f i-l ,3+1 + £ i+l,,j-l f i-l,. j-l 

dt ^ 2Ax 2Ax 

(2.15) 

g i+l,j+l ~ g i+l , j-l , g i-l,j+l ~ g i-l , j-l = 

2Ay 2 Ay 


We shall see that the stability properties of (2.15) are even better than 

(2.14) and allow an optimal time step. Now, however the even and odd points 
are decoupled. In the next section we discuss the accuracy of these schemes. 
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III. ACCRUACY - ONE DIMENSION 

We now consider the accuracy of schemes for nonuniform meshes. To 
simplify matters we first consider one-dimensional schemes. We consider the 
equation 

w + f =0. (3.1) 

L X 

We shall also consider the change of variables £ = £(x) and then (3.1) 

becomes 

f E 

w + — = 0 (3.2a) 

or in conservation form 

(x^ w) t + f ? = 0. (3.2b) 

£(x) will be chosen so that the nonuniform mesh in x becomes a uniform mesh 
in 5 . We will usually scale £ so that A£ = 1. Frequently £(x) is never 
constructed but the variable £ determines the location of x^ . In one 

dimension £ has the property that £(x^ ) = j at all the nodes, j = 0,«»«,N. 
To construct £(x) is now a standard problem in interpolation theory. 
However, we must have that £ 0 and would prefer that £(x) be 

monotone, i.e., if x^. x <_ x^ + ^ then j _< £(x) _< j + 1. This guarantees 

that intervals get mapped onto corresponding intervals. In multidimensions we 

would like to map each cell onto a standard rectangle or cube. If the 

geometry is so distorted that this cannot be done then this approach is not 


viable . 

We 

stress that 

there is 

no 

need to 

actually construct 

5 and hence 

there 

is no 

need for 

5 

to 

be 

defined 

globally. 

Since 

all accuracy 

requirements 

are local, 

it 

is 

sufficient 

that £(x) 

can be 

constructed 


locally at each mesh cell but with sufficient smoothness. 




Figure 2 

Before discussing accuracy we first mention how nonuniform meshes may be 
constructed. One difficulty with the analysis is that the theory of accuracy 
is an asymptotic theory as the mesh is refined. In practice one usually deals 
with a finite grid which is not extremely fine, especially in three 
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dimensions. Multidimensional grid generation about general configurations is 
difficult [16, 17], Not only is the mesh nonuniform but singularities can 

appear at which quadrilaterals collapse into triangles. One advantage of a 
finite volume approach is that the formulas are still valid in this limit. 

However, it is not clear what happens to the accuracy of the scheme in this 
degenerate limit. In Figure 2 we show a two-dimensional section of a three- 
dimensional grid showing a few cells near the leading edge. In the next plane 
triangular elements appear. This grid was generated using FL057 . 

Given a mesh there are two ways of constructing a finer mesh. One method 
is to subdivide each cell into smaller cells. This process usually leads to a 
smooth distribution but is difficult to implement using general 

quadrilaterals. An alternative is to grid the entire region anew using more 
nodes and ignoring the coarser mesh. This is equivalent to refining the mesh 
in the computational £ variables. With this method the mesh is less likely 
to be smooth in the physical space. 

Let Tj = hj + j/hj , with hj the local mesh spacing. We now introduce 

some notations that govern the rate of stretching. 

Definition ; 

(a) The grid stretching is quasi-uniform (or algebraic) if 


v. = 1 + 0(h p ), p > 0 


(3.3a) 


h 


max h . 


We note that the larger p is, the more smooth is the stretching 
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(b) The grid stretching is exponential if (a) is not valid but 


= 1 + 0(h). 


(3.3b) 


(c) The grid stretching rate is faster than exponential if neither (a) nor 
(b) are true. 

We note that (3.3a) can be reformulated as 


Vi - h j ’ 0<hP+1) - 


(3.3a") 


A typical mesh generator gives an a priori distribution of points in 
certain directions. Thus, in one space dimension one possibility is 


Xj = AK a j /N j = 0, N. 


(3.4a) 


Let hj = xj - Xj_j, we than find that 


(3.4b) 


a and A are chosen so that xq and x^ are fixed. Hence 


A = x Q * 0 


(3.4c) 


log(x N /x Q ) 

log(x) 


(3.4d) 
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Thus, for fixed k t 1 + 0(h) we have a fixed ratio between the areas of the 
cells and so the stretching is exponential. Furthermore, it is shown in [8, 
20] that exponential stretchings can accelerate the convergence to a steady 
state. On the other hand it is shown in [6] that nonuniform meshes may 
introduce false reflections from the gradients of the stretching. Exponential 
stretchings are common in many exterior calculations. When k = 1 + 0(h) 
then the stretching is algebraic. 

We shall only consider the case that the metric coefficients are evaluated 
numerically. Hoffman [7] has considered some examples where £ is a known 
function of x and all metric derivatives are calculated analytically. He 
has shown that if Ax is refined so that A£ is halved then quadratic 
convergence is retained. This is equivalent to having an algebraic grid 
stretching. 

We wish our formulas to be in conservation form. Again considering the 
integral form the equation we wish to solve is 

Jr + ^ / ?*ndS = 0 (3.5a) 

fi 

where w is the volume weighted average of w. We say that the scheme is in 
conservation form if we replace the integral by any consistent 

approximation. Hence, we replace (2.5) by 

dw. 

V j dt 1+ l \ F k = 0 (3 * 5b) 

J k 

where are weights that depend on the nonuniform spacing. In general 

will depend on several h j . We must however use the same summation formula in 


each cell 
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We first consider a finite difference mesh where both the dependent 
variable Wj and the geometric variable Xj are given at the nodes (see 
Figure 3a). 


>1 


w. 


w 


j+1 


x j-l 


X. 


x j+l 


he 


h. 

J 




K - 


h j+i 


-a 


Figure 3a 


Let hj = xj 


Xj_^. We consider the approximation to (3.1) 


Scheme I: 


dw . 


+ ct . f. , +6. f. + y • f • , i = 0. 

dt J j-1 j 'j j+1 


(3.6a) 


By Taylor series we find that this is second-order accurate if 


a . 
J 


. 1+1 


V h j + Vi> 


(3.6b) 




h . . - h . 
-1+1 3 

h. h. . 

1 J+1 


(3.6c) 


Y j 


Hence we have 


h J+ i <h j + 


(3.6d) 


dw. 

i 

dt 


h. 

J 


h. + h. . 
J J+1 


f.,, - f - 

(4^) + 

j+i 


h j+i 


h. + h. .. 
J J+1 


f. - f. , 

(-■ V — ) = °- 

j 


(3.6e) 


However, this scheme is not conservative 



-14- 


Using the same mesh as in Figure 3a we consider a finite volume scheme. 
In this case we replace (3.1) with 


d 

dt 


x j + l 

/ 


wdx + f 


j + 1 



(3.7) 


We next approximate the integral by W j^ x j+1 ~ X j-1^ which is only first- 
order accurate unless x^ is half way between x j-i and x j + l* then have 

Scheme II: 

dw. f . * - f . . 

+ JCL.o (3.8a) 


dt 


x. , - x. , 

j+i j-i 


or in terms of hj 



f j+l 


- f 


+ Vi 


o. 


(3.8b) 


We can also rewrite this scheme in a manner similar to that of (3.6e) as 


dw . 

3 

dt 


+ 


h. 

3 


Vl 


Vl 




') + 


h. 

J 


h. 

3 


+ V 


f. - f. , 

( ] h r ‘ ) - 0 . 

j 


(3.8c) 


Hence, Scheme II can be interpreted as a volume weighted average of the fluxes 

at j - I /2 and j + V 2 • The flux balance at j - V 2 an( l j + V 2 are correct 

to second-order. Nevertheless a Taylor series expansion verifies that this 

2 

scheme is first-order for a time-dependent problem unless h^ +1 - h_. = 0(h ), 
i.e., unless the mesh is quasi-uniform. In a steady state problem we would 
have to account for a source term or a second dimension so that the problem 
would be nontrivial. One can have an algorithm that is first-order accurate 
in space for a time-dependent problem but is second-order accurate in space 
for a steady state problem. This scheme is also conservative in the sense 


that 
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I ( 



w. 

J 


is constant except for boundary effects. Thus the conservation takes place 
if Wj is associated with part of the cell to the left and an appropriate 
part of the cell to the right of Wj . 

We next consider the mapping E, = £(x). A straightforward central 
differencing to (3.2a) leads to (3.8). We thus conclude that Scheme II is 
first-order (for nonuniform meshes) in x but is second-order in £. In 
other words, if w is quadratic in £ then (3.8) is exact but if w is 
quadratic in x then we have a discretization error. In theory we use a 
nonuniform mesh because it reflects the gradients of the solution. Hence, it 
would be natural to assume that w behaves quadratically in £ rather than 
in x. In practice the mesh is determined more by geometrical considerations 
than the behavior of the solution. Furthermore, for a system of equations not 
all the quantities grow at the same rate. When we are interested in global 
quantities, e.g., lift and drag, it is not clear if the approximations will be 
first- or second-order accurate. 

We next consider a modification of Scheme II. 


Scheme Ila: 


where 


dW j + f j+ V? f 1- V? , o 


£ j+v 2 ‘ “j £ j+l + (1 " “3 >f J 


a. = — s ^ — , h. = x. - x. , 

3 h 2 + h 2 , 3 3 3-1 

J 3 + 1 


(3.9a) 

(3.9b) 

(3.9c) 



(3.9d) 
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This scheme uses dependent variables at j - 1, j, j + 1 but geometrical data 
at j ~ 2, j - 1, j, j + 1. A symmetric variant is obtained if (3.9c) is used 
to define a j-i rather than a j . One can also obtain a similar formula that 
depends on the geometrical data at j - 2 through j + 2. This scheme is 
second-order in physical space as long as the mesh stretching is no worse than 
exponential. The scheme is conservative in that £ w^ is conserved. 

We now consider a finite volume scheme where the dependent variables are 
defined in the center of cells as shown in Figure 3b. This corresponds to the 
two-dimensional scheme described in the previous section and used in [8], 


W j-1 


w. 

J 


W j+1 


X j-3/2 


fc- 


Vi 


x j- v 2 


x . 


j + v 2 


he- 


j+i 


x j+3/2 


-$1 


We then replace (3.1) by 


Figure 3b 


V 1/2 


d? [ “ dX + f J+ Vi ' £ J- 1/2 ‘ ° 


X. 


- v 2 


(3.10) 


which as before is exact. Using the procedure described in Section II, e.g. 
(2.11a), we approximate (3.10) by 


dw. f + f . f . + f . , 

(v _ v i 3 . 3+1 3 _ 3 3~1 _ n 

3+V 2 j- V 2 ^ dt + 2 2 


or 
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Scheme III: 


dw . 

-J-J.+ 

dt 


£j+i £ j-i 

2(x j+V 2 - X j-V2 ) 


= 0 . 


(3.11) 


Assuming that f j is evaluated half way between Xj_ y and x^. + y^ , a Taylor 
series expansion verifies that this scheme is inconsistent unless 


h . , , - 2h . + h. , 

-Jii ^ 0(h), 

j 


(3.12) 


i.e., unless the mesh is quasi-uniform (algebraic) (3.3a). Note that (3.11) 
is second-order accurate only if the left-hand side of (3.12) is O(h^), i.e., 

n 

r j = 1 + 0(h) which requires an even smoother grid, i.e., (3.3a) with 

p = 2. Introducing a mapping E, = 5(x) we again verify that (3.11) is 

second-order in A£ . Hence, we conclude that for an irregular mesh (3.11) is 

inconsistent in x but is second-order in £. However, now the location of 

Wj is different in the two interpretations. Considering (3.11) in physical 

space Wj is evaluated half way between x j_]y and x ^ + y^ , i.e., 

x. = l/o (x . i, + x 1( ). However, in f space we have x. = x( i ) where 
J z J- 72 j+ V2 J VJ 

x = x(?) is the inverse map of £ = £(x), j - V 2 = S( x j_l/ 2 )> 

j + 1/2 = £(x. + ^). Hence, we assume that £(x) is invertible and that if 

j - V 2 I C < j + V2 then X j_ y < x(^) £ x^ + y . Since £ does not enter into 

the scheme we only need know that such a £(x) exists. The final scheme does 

not depend on the derivation of (3.11) as coming from a physical or 
computational space. However, the final interpretation of the numerical 


solution, e.g., graphics, does depend on how one defines x j • Furthermore, if 
the fluxes depend explicitly on the geometry, e.g., cylindrical coordinates, 
then the two interpretations of (3.11) will result in different schemes. 
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In order to remove the inconsistency of (3.11) for an exponential mesh we 
consider another approach previously used by Collela and Woodward to construct 
upwind differences [5], Let 


x 


W(x,t) = 


Then 


3 + V 2 j“ l k x ^_ 1 /, 


f w(x,t)dx. 

j- \ 


" ( V v 2 5 ' ", 


(3.13) 


(3.14a) 


is the volume average of w used in the finite volume scheme (3.10). We also 
have that 

W(x j _ i /2 ) =0 (3.14b) 


W(x. +V? ) = w. + w . 

J+3/ 2 J x j + v 2 " x j-v 2 J+1 


- x 4+l/ 


(3.14c) 


The local value of w(x) is then given by w(x) = We now fit W(x) with 

O X 

a quadratic, interpolating at x. , x. i. , x_ Q/ „ using (3.14). We then 

3 /2 3+ h J + d/2 

differentiate W(x) and evaluate at Xj + 1^ to find 


W. 


j + V 2 


h . W . , i + h . , . W . 
... 3 J + 1 J + 1 3 


h. + h.,, 
3 3 + 1 


(3.15) 


Then (3.9) becomes 
Scheme Ilia: 


dW . f (W. , i/„ ) — f (W . 


dt 


l A. 


\zlhL = o. 


S j + V 2 X 3 - V 2 


(3.16) 


Interpreting W^ as w at the midpoint between x j- 1/2 anc * X j+ V 2 
expanding in a Taylor series we find that (3.16) has a truncation error 


and 
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t = 0(h. +1 “ h j) + 0(h ). 


Hence, (3.15), (3.16) is second-order accurate for an algebraic mesh and 

first-order accurate when used with an exponential mesh. Thus fj + is 
simply given by linear interpolation in W accounting for the nonuniform 
mesh. We note that if we wish (3.10) to be fourth-order in space we do not 
wish to find f j + 1^ using cubic interpolation. This occurs since fj + i j 


- £ 


d£ 


. i. is only a second-order approximation to - — . Hence, to get fourth- 
J“ 12 dX 

directly using x j-2 through x j+2* This 


order accuracy we must find 


i! 

3x 


has been previously observed by Zalesak [23]. We also note that the scheme is 
conservative since (3.10) is derived from the integral form. In particular 


| -f h j "j 

is constant in time excluding boundary effects. The use of (3.16) gives 
first-order accuracy if r^ is almost constant, i.e., the mesh is 
exponential, and gives second-order accuracy if the mesh is algebraic, i.e., 
rj = 1 + 0(h). 

We wish to stress that the difficulties with the cell centered quantities 
arises because of the finite volume approach. If we define the dependent 
variables in the center of the cell but use a finite difference approach then 
we can achieve higher accuracy as was done for the case of nodal variables in 
Schemes II and Ila. In particular we consider the following generalization of 


Scheme Ila 
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Expanding 


« -77 'h Vi + (1 - “j ,£ j] - h-i h - (1 - 


<3 


in a Taylor Series we find that 


Hence, 


'! = 77 ( l /2 + h j + i) + U - + hj-i > ] f ' 

+ ^ [Oj thj + h. +1 ) 2 - (1 - a^Xh. + h.. 1 ) 2 ]f" + 0(h 3 )}. 
Q is a first-order approximation to f' if 


“j - l/2 h (h 3 


+ h j+i ) + 


(1 - Oj.pO.j + 


Vi>] 


+ 0(hP), p = 2. 


(3 


Furthermore, Q is a second-order approximation if (3.19) holds with p 
and in addition 


cXj(tu + h j +1 ) 2 = (1 - Oj.jKhj + h j_^) 2 + °(h 3 ) 


(3 


or letting r^ = h j +1 / h j 


a.(l + r.) 2 = (1 - a. ,)(l + 2 + 0(h). (3. 

J 3 J r j_l 

For a finite volume scheme k. = h. and so (3.19) determines n and we 

J J “j 

Scheme Ilia which is first-order. For a finite difference scheme we 
choose both a. and <., We thus have 


.17) 

.18) 

.19) 
= 3 
i.20) 

20a) 

have 

can 
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Scheme Illb: 



f j +1 /7 f .i -V? = 

K . 

J 


(3.21a) 


= a. + (1 - a.)f . 

‘ J J 


a . = 


j+ V2 " “j Vi 
(h j + h j-i )2 


j (h. + h. ) 2 + (h. + h ) 2 ’ j j+1/2 2 1/2 

J J-l J J+l 


(3.21b) 


(3.21c) 


(h + h ) (h + h ) 

k — a — “ J + ( 1 — ct ) — * * 

j a j 2 u a j-i ; 2 


(3.21d) 


This scheme is now second-order in physical space as long as the mesh spacing 
is no worse than exponential. However, viewed as a finite volume scheme 
(3.21a) implies that 

V v 2 

J wdx = k. w. (3.22) 

V 1/2 J J 


which is not within the spirit of a finite volume approach especially for 
multidiraensions . 

We now summarize the conclusions of this section. In general both central 
finite difference formulas and finite volume formulas will be only first-order 
accurate in the physical variables when nonuniform meshes are used unless the 
weighted formulas (3.6) are used. When the dependent variables are defined in 
the middle of the cell then the formulas are not even consistent unless the 
mesh is quasi-uniform. In most cases the scheme is second-order measured in 
terms of some computational variables £(x). If this variable £ is 
carefully chosen to match properties of all components of the solution, then 
this second-order accuracy is meaningful. However, if the stretching is 
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utilized simply to put the outer boundary far away then the stretching can 
reduce the order of accuracy. Numerical experiments verify that the 
computational solution can be very inaccurate although formally the solution 
is second-order in £. One could possibly post-process the solution, but this 
would be very difficult to do on general multidimensional meshes. 

In this section we have only discussed seraidiscrete formulations where the 
time discretization is not included. It is easy to see that the extensions to 
fully discrete algorithms using Runge-Kutta type formulas in time, e.g., [8, 
20] or else implicit formulas, e.g., [2, 3, 11], are straightforward. When a 
MacCormack type formula is used then one can again verify that one loses 
accuracy on a nonuniform grid. In this case to achieve second-order accuracy 
in x one should let h. = x. - x. , and then 


and 


n 


2h. 


w. = w. - 7- — ■ 

J J h. + h. 


At. 


j+1 


ft* - f n 

f J + 1 J 

Vi 


(3.23) 


n+1 
w . 

J 



+ w. 
J 


h. 

J 


2h 


3 + 1 


+ h j + l 




IV. ACCURACY - TWO DIMENSIONS 

In one dimension we have seen that when the geometrical data and the 
variables are defined at the same point then a simple averaging of the fluxes 
gives first-order accuracy in x and second-order accuracy in the 
computational space variable, £. When using nonuniform meshes the accuracy in 
the physical space can be improved by using the weighted formula (3.6). When 
the dependent variable is defined in the middle of a cell then simple 
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averaging of the fluxes can give an inconsistent scheme for an exponentially 
stretched mesh. Instead a weighted average (3.14), (3.15) should be used. 

In extending these results to two space dimensions we encounter several 
new difficulties. First a x difference at given point i can be given by 
evaluating it at (i, j) or averaging over several j lines, e.g., (2.13), 
(2.14), (2.15). In addition weighting formulas, e.g., (3.14), (3.15) must now 
be replaced by two-dimensional interpolations. By considering cells with 
equal volume but different shapes it is easy to see that volume weighting can 
decrease the accuracy of the algorithm. It is necessary to include more 
information than just the volume of neighboring cells. 

In one space dimension we defined a smooth (or quasi-uniform) grid as one 


for which r. 

J 



1 + 0(h). In two space dimensions we need to consider 


two-dimensional effects and not just demand that the volumes vary smoothly. 
Consider the situation shown in Figure 4. All the cells have the same area 
but it is obvious that one should not define the flux at G as the average of 
those at the cell centers E and F. Rather the flux at G must be found by 
a two-dimensional interpolation. 



Figure 4 
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If we wish the standard averaging to be higher order for sufficiently smooth 
grids, then our notion of smoothness must include changes in angles between 
zones as well as changes in volumes between zones. We thus define a 
multidimensional grid as being quasi-uniform if the change between cells of 
both x and y along each side is small. For example in Figure 4 the grid 
is quasi-uniform only if 


= 1 + 0(h) 


= 1 + 0(h), 


(4.1) 


and similarly for the other three sides. It is easy to see that (4.1) holds 
only when the relative change in both the volume and in the angles is 
1 + 0(h). 

We first consider the generalization of Scheme II where the dependent 
variables and the metrics are defined at the nodes (i, j). We perform the 
finite volume integration over the four cells shown in Figure 5. 



A B C 


Figure 5 
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We then have 

4 - // wdv + / (fdy - gdx) = 0, (4.2) 

« r i +r 2 +r 3 +r 4 

where the line integral is performed in the clockwise direction. We 

approximate the area integral by w. . »V, where V is the area of the four 

1 » J 

cells. For the line integrals we can use either the midpoint or trapezoidal 
rules. The trapezoidal rule is second-order accurate while the midpoint rule 
is only first-order since (i + 1, j) is not necessarily at the midpoint of 
T 2 * Since the final formula will only be first-order accurate it is not clear 
if this makes a difference. Using the midpoint we approximate (4.2) by 


dt + f D Ay EC g D AX EC ~ f F Ay EC + g F 


" f H Ay GH + g H AX GA + f B Ay CA ' g B AX CA = ° 


where V(fi) - l / 2 (Ay EA Ax CG - Ax EA Ay CG ) and A x CA = x c - x A etc. 

We see from Figure 5 that H is not necessarily near to being the 
midpoint of T^. A better approximation to an integral along is to use 

the trapezoidal rule on each section. Thus, for example, 

f r + f H f A + f H 

J fdy = (Ay) HG + 2 - ' (Ay) AH* (4 * 4) 

Using a uniform grid (4.3) would reduce to (2.13) while (4.4) would reduce to 
(2.14). An intermediate possibility is to use the trapezoidal rule on all of 
r 4 . We then would approximate 
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f r + f a 

J fd y s 2 Ay AG* (4,5) 

For a uniform mesh (4.5) would reduce to (2.15). We refer the reader to the 
appendix for the full two-dimensional formulas. 

In this discussion we have not considered the effects caused by the 
differencing of metrics. In the far field exterior to a body the flow is 
almost uniform but the mesh is frequently highly stretched. It is desirable 
that the derivatives be small independent of the mesh stretching. Steger [14] 
discusses various techniques to accomplish this. It is obvious that finite 
volume schemes automatically satisfy this grid conservation property, see, 
e.g., (A2), (A4). In addition Steger also presents numerous examples that 

demonstrate the effect of nonuniform or highly skewed grids on the accuracy of 
the solution. For some schemes highly stretched meshes can also affect the 
rate of convergence to a steady state. 


V. TIME INTEGRATION 

All the formulas discussed until this point have used a semidiscrete 
formulation and have concentrated on the space discretization. One possible 
way of integrating these equations in time is to use an implicit method with 
A.D.I. splitting for multidimensional problems (see [2], [3]). In the rest of 
this section we shall describe and analyze the use of Runge-Kutta type 
formulas [8]. 

Let Lu represent the space discretization of a method. Then a K stage 
Runge-Kutta method has the form 
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(1) _ n . 
u v ' = u + a, 


AtLu 


(k) _ n (k-1 ) . „ 

u - u + AtLu k = 1,2, ••• > K 


(5.1) 


where u n+1 = u^ . are given parameters with = 1 for at 

least first-order accuracy in time. When L has constant coefficients we can 
Fourier transform (5.1). Let Q ( 0 , ip) be the Fourier transform of the space 
discretization operator L and let u be the Fourier transform of u. We 
then find that (5.1) becomes 

u n+1 = Gu n (5.2a) 

where 

G = P R (QAt) (5.2b) 


and P„ is a polynomial 


K 


P K (z) = l 0. z. 


j=0 


J J 


(5.2c) 


where 


0 O 1 > a K> 6 j a K a K-r*” ,a K-j+r 


(5. 2d) 


(5.2c) generates a stability region, i.e., a domain in the complex plane 
where |P^,(z)| <. 1. We shall now consider the hyperbolic (inviscid) case. 
All the formulas introduced have central differences and so Q lies on the 
imaginary axis. In this case (5.1) is stable if 


At <4 

A 


(5.3) 


where K depends on the sequence a,, •••,(*„ and 

1 K. 
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X = max p(Q(0,iji)) (5.4) 

— it <9 , 

A A 

where p(Q) is the spectral radius of Q. Hence, given a sequence 

A 

the allowable time step depends on the eigenvalues of Q, i.e., on properties 
of the scheme. 

Specifically we consider the Euler equations in general curvilinear 
coordinates. This system can be written as 


w + Aw + Bw = 0 (5.5) 

t x y 


where (x, y) are general coordinates w 


(p, pu, pv, E) 


t 


and 



Y 

y 

q - (y-2)Y y u 

Y v + (y-l)X u 
y 1 x 

Y y h - (y-l)qu 
-Y 

x 

r + (y-2)Y u 
x 

-Y v - (y-l)X u 

X 1 X 

-Y h - (y-l)ru 
x 


-X 

y 

-Xu- (y-l)Y v 

y y 

q + (y-i)x x V 

-X y h - (y-l)qv 
X 

x 

X u + (y-1)Y v 
x x 

r - (y-2)X x v 
X x h - (y-l)rv 
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Here X, Y are the Cartesian coordinates and (x, y) are the generalized 


coordinates. Also 


q = Y u - X v 

y y 


r = X v - Y u 

X X 


(5.7) 




c, 2 2 

u _ 2 t u + v 

h i- + 5 . 

y - 1 2 


We first consider the case that 


„ £ i+l,j 


g - g i,j + l ~ g i,j-l 
’ 8 y 2 Ay 


(5.8) 


This scheme (2.11a) is used in FL052, [8] and reduces to (2.13) for a uniform 
mesh. In this case the stability criterion becomes [20] 


(5.9) 


|q| + |r|+/x 2 +Y 2 +X 2 +Y 2 + 2|x X + Y Y Ic 

» * l l I v v \r • \r v l 


x x y y 


x y x y 1 


We next consider the case where the flux integrals are approximated by the 
trapezoidal rule. This scheme is given in (2.11b) and reduces to (2.14) for a 
uniform mesh. The schemes given in the appendix, which account for 
nonuniformities in the mesh while averaging, are also included in this case. 


We thus consider 
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f j+l,j+l f i-l,j+l + 2 ^ f i+1 , j f 1-1 , j ^ + f i+l,j-l f i-1 , j~l 

8 Ax 


(5.10) 


S-.-L1 S-i 


g y = 


i+l,j+l 6 i+l , j 


,j-l + 2(S i,j+l " + g i-l,j+l 


8Ay 


For this case the stability limit is given by [19] 


At < 


K 


I w I + /a 2 ■ ’ 2 


/a + b c 


(5.11) 


with 


Y + Y + max(Y , Y ) 
x y 


a --5 y - 


X + X + max(X , X ) 
v \t x x 9 y ' 


b = -£ * 


(5.12) 


„ _ q + r + max(q, r) 
W 2 


We note that this is a larger time step than allowed by (5.9) but at the 
expense of a more complicated algorithm. In FL052 the calculation of the 
fluxes and the artificial viscosity accounts for most of the running time. 
Hence, complicating the averaging procedure for the fluxes does not require 
much extra computer time. 

Finally we consider the averaging 

f * (f i+l,J+l ~ f i-l,3+l ) + (£ i+l,j-l ~ 

x 4 Ax 

s (g i+l,j + l ~ gj+lj-^ + (g i-l,j + l ~ ^-l.j-^ 

K y 4 Ay 


(5.13) 
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In this case we can give two criteria that are sufficient for stability (see 
also [1]). One is that 

At < K max [^j, ^|y] (5.14) 

or equivalently 


At « ——— . (5.13) 

maxMql + c, I r I + cl 

111 y y 1 1 x x J 

This is the same condition one would obtain using time splitting. A second 
possible criterion is that 

At < K (5.16) 

max( | q | , | r | ) + /max(X^ + X^, Y^ + Y^)c 

Both these criteria provide sufficient conditions for stability. In both 
cases these allow larger time steps than (5.11) which in turn allow larger 
time steps than (5.9). 


VI. STABILITY 

We next consider the effect of a nonuniform grid on the stability of an 
explicit scheme. To be specific we will compare the Schemes III and Ilia for 
the one-dimensional problem. For an equation with constant coefficients 


Scheme III (3.10) becomes 


w^ + aw = 0 
t x 



Qw = 0 


(6.1a) 
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with 


Qw = 


a( Vi - “1-1 > 




(6.1b) 


Formally replacing Wj by e 1 ^ (i.e., forming the pseudodifference 

operator) we find that 


Q = 


ai 


X j+V 2 " X j-V 2 


• sin £. 


( 6 . 2 ) 


We thus see that the symbol is purely imaginary. Using a Runge-Kutta type 
scheme [8] (see also (5 . l)-(5.4) ) , this leads to a stability limit of the form 


h j = x j +1 /2 x .i~ V?. 


At. 

J 


At . 
J 


< K/a 


(6.3) 


where K is given by the intersection of the stability curve of the Runge- 
Kutta scheme with the imaginary axis. 

We next consider the Scheme Ilia (3.14)— (3. 15) with constant coefficients, 


where 


dw. 



w 


0 


w. 

J 


( h . w . , , + h . w . 

-J i±i J +1 3 

h. + h 
J J+1 


h . . w . + h . 
J- 1 J 3 
h. . + h. 
J-l J 



(6.4a) 


(6.4b) 


In this case 


Q a =E7 


h j (h j+i ' 


(h. + h j+1 )(h._ 1 + h.) 


(1 - COS ?) 


+ ill + 


h 1 ~ h 1-l Vl 

(h. + + h.) 


Ism 


(6.5) 
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Note: Q has a nonzero real part. Q is dissipative if and only if 


a(h J+1 " VP > °* 


( 6 . 6 ) 


In other words, if the mesh expands in the direction of the wave then the wave 
is dissipated. If the mesh contracts in the direction of the wave motion then 
the wave is amplified. Alternatively, consider an 0 mesh where the spacing 
between nodes increases as the mesh goes to the outer boundary. Then outgoing 
waves are dissipated while incoming waves are amplified. Thus absorbing 
boundary conditions are important. We point out that the above statement is 
not completely accurate. The stability region of many Runge-Kutta schemes 
include a portion of the positive complex plane. Hence, even though the real 

A 

part of -Q is positive, nevertheless, the total Runge-Kutta scheme can be 
dissipative. Usually this intrusion of the stability region into the right 
half plane is largest near the imaginary stability limit. Hence, the 
amplification of incoming waves is least pernicious when the time step is near 
the stability limit. This is an additional motivation to use a local time 
step. 


The time step restriction for (6.5) is much more difficult to analyze than 

A 

for (6.2) since Q contains real and imaginary parts that depend on the 

Fourier variable £. As before we let r. = h... ,, , then 

J 3 +1/h j 


Q a = r~ 
a h . 

J 


r. r. . - 1 
3 3-1 


(1 + rTKl + rT_j) 


(1 - COS £) 


+ i 1 + 


r. . - r. 
3-1 3 


(1 + r )(1 + r x ) 


sin £ 


(6.7) 
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Hence, if the mesh is algebraic then Q = Q + 0(h) and so the stability 

3 . 

limit is given by (6.3). When the mesh ratio is exponential then the 

A 

imaginary part of Q (i.e., the wave speed) is changed by 0(1) terms. 
Hence, for an exponentially stretched mesh the stability properties of the 
scheme are different than that for a uniform mesh. 

In this entire analysis we have assumed that the coarsest mesh is still 
sufficiently fine to resolve the waves. Hence, stability analysis,, which is 
an asymptotic theory, is still valid. In many cases the outer meshes are 
large compared with the relevant wave lengths. In this case other phenomena, 
e.g., wave reflections, may also be present, see [6], 

We next consider the effect of using a local time step for the fully 
discrete algorithm. (6.1) can then be approximated by 


,n+l „n 

V . 

J 


aAt . 


WV‘ - wV + — , 1 7 (w n , , - w n ) = 0 

J W J+1 J_1 


( 6 . 8 ) 


choosing At^ = K(x^ + ^ - x^_ y ) as the local time step (6.8) becomes 


W n+1 - W n + ^ (W?. , - w n , ) = 0 
J J 2 ^ j+1 j-l' 


(6.9) 


which is stable when used in a Runge-Kutta type mode for the appropriate K. 

However, now all effects of the nonuniform mesh disappear. Hence, in one 
dimension we have recovered the standard analysis. Thus using a local time 
step removes the dissipation or growth previously discussed. It also removes 
the false reflection found by Giles and Thompson [6]. Using the finite volume 
scheme with linear interpolation, (6.4), the use of a local time step cannot 
remove the real part of the eigenvalues and so we cannot recover the 
properties of a uniform grid. 
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VII. COMPUTATIONAL RESULTS 

In order to show the value of accounting for the nonuniform mesh we 
present one sample case. We consider the Runge-Kutta scheme for the Euler 
equations, FL052T described in [8], This is a finite volume method wit^i the 
unknowns given at the cell centers. As a specific case we consider the 'two- 
dimensional flow about a NACA0012. The inflow conditions are M = 0.3, * 

CO 

a = 10°. Since this is a subsonic flow the numerical solution can be compared 
with the solution to a potential code with a very fine mesh. In Table I we 
present the lift and drag for a coarse 64 x 16 0 mesh. The original code is 
based on a two-dimensional version of (3.9)— (3.10) while the new version is 
based on (3.14)— (3. 15) . Though this correction only accounts for one- 
dimensional effects we note a marked improvement in the results. 


Table I 


Code 

C L 

C D 

Potential code 



with fine mesh 

1.274 

0.0000 

Original FL052T 



64 x 16 0 mesh 

1.227 

0.0087 

Improved version 

1.264 

0.0025 


Comparison of new and original code for flow about NACA0012 with = 0.3, 
a = 10°. 
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VIII. CONCLUSION 

In (3.3) we defined algebraic and exponential rates of grid stretching. 
With an algebraic stretching, the grid metrics are sufficiently smooth that 
all second-order techniques retain their formal accuracy. When exponential or 
faster stretching is used then the accuracy of a central difference scheme or 
finite volume method will deteriorate to first-order accuracy unless special 
weighting formulas are used, e.g., (3.6). When the dependent variables are 
defined at the center of the cells, then the finite volume approach can yield 
an inconsistent scheme unless the weighting (3.14) is used. With this 
weighting the scheme is second-order for an algebraic stretching and is first- 
order when the stretching is faster than algebraic. 

In all cases the algorithm is second-order measured in terras of a 
computational variable § rather than in terms of the physical coordinate 
x. If this 5 corresponds to the behavior of all components of the solution, 
then this measure of accuracy is reasonable. It is clear that £ should 
always be chosen so that the solution is smooth in the £ variable. However, 
in many cases £ is chosen simply to put the outer boundary sufficiently far 
away. Thus, for example, if the solution has a boundary layer then one would 
want the grid to be almost uniform within the boundary layer. One could use 
an exponential mesh in the outer region only if the solution decays 
exponentially in the far field. Furthermore, if £ = £(x) is not given 
analytically but rather xj is given then there may not exist a sufficiently 
smooth and monotone £(x). 

In multidimensions these difficulties are compounded. First the grid is 
quasi-uniform only if both the volumes and all angles vary sufficiently 
slowly. Difficulties in constructing multidimensional grids leads to 
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nonuni form grids that are due to geometrical effects that have nothing to do 
with the behavior of the solution, e.g, Figure 2. Hence, in this case, it 
becomes more important to introduce weights that compensate for the nonuniform 
grid. Unfortunately the weighting formulas must now account for the 
multidimensional geometrical behavior. One possibility is to use bilinear 
interpolation within a cell as is done in finite element codes, see formulas 
in the appendix. 

In several dimensions the line integrals of the fluxes can be approximated 
by either a midpoint or a trapezoidal rule. For a uniform mesh with cell 
centered variables the midpoint rule yields a simpler formula and one that 
gives no difficulty for the tangential flux near a boundary. However, when 
the mesh is distorted the midpoint rule may be more inaccurate due to 
difficulties in locating the actual midpoint. Furthermore, if a two- 
dimensional interpolation is used then the simplicity of the midpoint rule is 
lost. The trapezoidal rule requires more work but allows a larger time 
step. In addition it is now easier to use either volume weighting or bilinear 
interpolation to calculate the value of the fluxes at the nodes. Weighting 
formulas similar to (3.14) can now be used; see the appendix for full two- 
dimensional formulas. Alternatively one can use a finite difference or finite 
volume scheme where all variables are defined at the nodes. This guarantees 
first-order accuracy in x even when the mesh is exponential. For two- 
dimensional problems, grids can be constructed which are fairly quasi- 
uniform. In this case both approaches give similar accuracy, see, e.g., 
[17]. For three-dimensional problems about complex bodies it is much more 
difficult to construct a quasi-uniform grid. In this case some sort of 
patching of grids may be necessary. 
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The use of a trapezoidal integration rule instead of the midpoint rule 
allows a larger time step but at the expense of a slightly more complicated 
formula and additional difficulties near the boundaries. The use of weighted 
formulas can change the dissipation and stability properties of the scheme. 

We also show that nonuniform meshes may change the dissipation properties 
of the time integration scheme. These dissipation terms are further changed 
if a local time step is used. 
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APPENDIX 

In this appendix we give the details for two separate two-dimensional 
schemes. The extension to three dimensions is immediate. We only consider 
the semidiscrete problem with the time integration to be explicit Runge-Kutta 
or else some implicit scheme. 

In the first algorithm we assume that both the coordinates and variables 
are given at the nodes (see Figure Al). We use a finite volume approach and 
evaluate the line integrals using a trapezoidal rule. 


G F e 



K 

(i»j) 



ABC 
Figure Al 


Consider the equation 

w + f + g =0. (Al) 

t x y 

We convert this to an integral relationship as in (2.6). Let f^ denote f 

at the point A and let Ax,,,. = x„ - x etc. for all the possibilities. 

D i 1 D H 

Note that the order of the indices is important. We then have 
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Scheme A1 : 


V ij dt + 1/2 t f A Ay BH + f B Ay CA + f C Ay DB + f D Ay EC 
+ f E Ay FD + f F Ay GE + f G Ay HF + f H Ay AG^ 

“ t S A AX BH + g B AX CA + g C AX DB + g D AX EC 

+ g E Ax FD + g F Ax GE + g G Ax HF + g H Ax AG^ = ° 


(A2) 


where 

V i,j =1/2 ^ Ay EA Ax CG " Ax EA Ay CG I * 


For a uniform Cartesian mesh this reduces to 


dw i,,j , f C ~ f A + 2(£ P ~ f H } + f E ~ f G 
dt 4 Ax 


. g G " g A + 2(g F ‘ g B ) + g E " g C _ 

-f- i u < 


(A3) 


4Ay 


The second scheme that we consider has the geometrical quantities at the 
nodes but the dependent variables are given at the cell centers (see Figure 
A2) this extends the scheme used in FL052 [8], However, we use the 
trapezoidal rule to evaluate the integrals rather than the midpoint rule used 
in [8]. 
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H r 

1 

1 

1 

1 C 

1 

1 

1 

1 

1 

1 

(i» j) ! 

l 

1 

1 

1 

E 

F 


A B 


Figure A2 


The results in 
Scheme A2: 


dw . 

V i,j dt ,J + [ (f C " f A^ Ay DB + B “ f D^ Ay CA^ 

(A4) 

- V 2 [ ( § c " S a )Ax db + (g g - § d )Ax ca ] - 0. 

In order to evaluate the ‘fluxes at the nodes we could use volume weighted 
averages. An alternative is to use bilinear interpolation using the cell 
centers. We first must locate the cell centers. Using simple averaging we 
find that 

x i,j =1/ * (x a + *B + X C + V 

(A5) 

y iJ =1/ * (y A + y B + y C + V* 




-42- 


As an example of using bilinear interpolation we give the formulas for 
calculating f at the point C given the values of f at the points E, F, 
G, and H and given the locations of all the nodes and cell centers. Let 


G 1 (f i+1 , j f i,j )AX GE (f i+l,j+l f i,j )AX FE 


G = (f .- f. . )Ax^ - (f . - f. . )Ax 
2 i+l, J i,3 HE i, J+l i,3 1 


FE 


a Ay pE Ax G£ Ay GE Ax pE 


(A6) 


3 = (x.,. . y.,, . - x. . y. . )Ax^,, - (x. Ll ... y.,. ... - x. . y. .)Ax„„ 
l+l, j 7 i+l,j i,j i,J Gh 1+1, 3+1 1+1, 3+1 i,3 i,3 FE 


y ■ 4y FE 4x HE - 4y KE iX FE 


6 * " v i,J y l,j )4x HE - (x i, j+ l y l,J« - x i,j y l,j> 4x FE . 


We then define 


6 G^ - 3G2 
a5 - 3y * 


aG 2 ~ yG 1 
a6 - gy 


_ f i+l ,i - f i„i - C4y FE - D(x i+l,j y i+l,j - X l„i y 1.1> 


Ax. 


FE 


(A7) 


A = f . . - Bx . . - Cy . . - Dx . . y . . . 
i,J i,J i,J i,3 1,3 


Finally we have 

f c = A + Bx c + Cy c + Dx c 


(A8) 
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For a uniform mesh Ax = Ay this gives 


(A9) 


For both Scheme A1 and Scheme A2 the stability condition is given by (3.11). 
This allows a larger time step than the stability criterion (5.9) for the 
scheme used in FL052. 
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